fert.data <- read_csv(here::here("fert-data.csv"))
## Parsed with column specification:
## cols(
## .default = col_character(),
## pH.experim = col_double(),
## Perc.Fertilization = col_double(),
## Insemination.mins = col_double(),
## Fert.success.mins = col_double(),
## Sperm.pre.exp.time = col_double(),
## egg.pre.exp.time = col_double(),
## pH.control = col_double(),
## pH.delta = col_double(),
## Sperm.per.uL = col_double(),
## Sperm.per.mL = col_number(),
## n.females = col_double(),
## n.males = col_double(),
## Yeear = col_double()
## )
## See spec(...) for full column specifications.
fert.data.2 <-
fert.data %>%
dplyr::select(pH.experim, Perc.Fertilization,
Insemination.mins, Fert.success.mins, Sperm.pre.exp.time,
egg.pre.exp.time, pH.delta, Sperm.per.mL, sperm.egg, n.females, n.males) %>%
mutate_if(is.factor, as.numeric) %>%
mutate_if(is.character, as.numeric)
## Warning: NAs introduced by coercion
fert.data.3 <-
fert.data %>%
mutate_at(c("Phylum", "Common name", "Brooders/Spawniers", "Family", "Taxa", "Species") , as.factor)
fert.data.3$pH.group <- cut(fert.data.3$pH.experim,c(6,7.3,7.5,7.8, 8.2))
fert.data.3$Phylum <- factor(fert.data.3$Phylum, levels = c("Echinodermata", "Mollusca","Cnidaria","Crustacean"))
str(fert.data.3)
## Classes 'spec_tbl_df', 'tbl_df', 'tbl' and 'data.frame': 465 obs. of 26 variables:
## $ Phylum : Factor w/ 4 levels "Echinodermata",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ Common name : Factor w/ 31 levels "amphipod","Antarctic bivalve",..: 10 10 10 10 19 19 23 23 21 21 ...
## $ Brooders/Spawniers: Factor w/ 3 levels "Brooder","Spawner",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ Family : Factor w/ 27 levels "Acartia","Acroporidae",..: 2 2 2 2 18 18 18 18 18 18 ...
## $ Taxa : Factor w/ 13 levels "abalone","amphipod",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ Species : Factor w/ 42 levels "Acartia erythraea",..: 4 4 4 4 16 16 36 36 33 33 ...
## $ pH.experim : num 8.15 7.85 8.15 7.85 8.15 7.85 8.15 7.85 8.03 NA ...
## $ Perc.Fertilization: num 99 90 42 48 87 86 97 98 57 43 ...
## $ Insemination.mins : num 180 180 180 180 180 180 180 180 420 420 ...
## $ Fert.success.mins : num 180 180 180 180 180 180 180 180 420 420 ...
## $ Sperm.pre.exp.time: num 0 0 0 0 0 0 0 0 NA NA ...
## $ egg.pre.exp.time : num 0 0 0 0 0 0 0 0 NA NA ...
## $ pH.control : num 8.15 8.15 8.15 8.15 8.15 8.15 8.15 8.15 8.03 8.03 ...
## $ pH.delta : num 0 -0.3 0 -0.3 0 -0.3 0 -0.3 0 NA ...
## $ Sperm.per.uL : num NA NA NA NA NA NA NA NA NA NA ...
## $ Sperm.per.mL : num 1e+05 1e+05 1e+05 1e+05 1e+05 1e+05 1e+05 1e+05 3e+06 3e+06 ...
## $ sperm.egg : chr "50000" "50000" "50000" "50000" ...
## $ n.females : num 1 1 1 1 1 1 1 1 3 NA ...
## $ n.males : num 1 1 1 1 1 1 1 1 3 NA ...
## $ Other : chr "used engauge digitizer" "used engauge digitizer" "used engauge digitizer" "used engauge digitizer" ...
## $ Data source : chr NA NA NA NA ...
## $ Notes/other treat : chr "temp = 27" "temp = 27" "temp = 27" "temp = 27" ...
## $ DOI : chr NA NA NA NA ...
## $ Author : chr "M. Schutter, Y. Nozawa, H. Kurihara, 2009" "M. Schutter, Y. Nozawa, H. Kurihara, 2009" "M. Schutter, Y. Nozawa, H. Kurihara, 2009" "M. Schutter, Y. Nozawa, H. Kurihara, 2009" ...
## $ Yeear : num 2015 2015 2015 2015 2015 ...
## $ pH.group : Factor w/ 4 levels "(6,7.3]","(7.3,7.5]",..: 4 4 4 4 4 4 4 4 4 NA ...
fert.data.4 <- fert.data.3[,c(1,5,7:9,11,12,16:19)] %>%
mutate_if(is.factor, as.numeric) %>% #convert factors to numeric factors
mutate_if(is.character, as.numeric) #convert character column to numeric
## Warning: NAs introduced by coercion
Gowers allows for missing data and multiple data types
dist.gower <- vegdist(fert.data.4, "gower")
Its usage is: cmdscale(d, k, eig = FALSE, add = FALSE) where: • d is a dissimilarity object (generated by dist or vegdist) • k is the number of principal components (PC) that should be extracted from the distance matrix (max number = min(col, rows)-1) • eig, logical. If TRUE eigenvalues for each PC are retuned. Default: FALSE. • add, logical. If TRUE a constant is added to each value in the dissimilarity matrix so that the resulting eigenvalues are non-negative. Default: FALSE.
The principal scores are contained in spe.pcoa$points and the eigenvalues are contained in spe.pcoa$eig.
spe.pcoa <- cmdscale(dist.gower, eig=T, add=T, k=2)
head(spe.pcoa$points, n=15)
## [,1] [,2]
## [1,] -0.5743930 -1.0530913
## [2,] -0.6098390 -0.8683141
## [3,] -0.7893031 -0.2945215
## [4,] -0.7845755 -0.2677329
## [5,] -0.6166947 -0.9347762
## [6,] -0.6298105 -0.8185384
## [7,] -0.5788613 -1.0403509
## [8,] -0.5844064 -0.9437828
## [9,] -1.0266021 -0.3297779
## [10,] -1.6639959 0.2962041
## [11,] -1.0892730 0.0416065
## [12,] -1.8912539 0.5178418
## [13,] -1.9568934 0.9443114
## [14,] -1.9565642 0.9317737
## [15,] -0.8084025 -0.9088631
head(spe.pcoa$eig, n=15)
## [1] 327.84640 193.43571 136.24688 87.00823 75.00325 68.77631 67.36747
## [8] 61.84909 58.84786 53.01254 51.42648 48.04094 44.58181 43.14809
## [15] 40.83116
hist(spe.pcoa$eig/sum(spe.pcoa$eig)*100)
plot(spe.pcoa$eig[1:100]/sum(spe.pcoa$eig)*100,type="b",lwd=2,col="blue",xlab= "Principal Component from PCoA", ylab="% variation explained", main="% variation explained by PCoA (blue) vs. random expectation (red)")
lines(bstick(100)*100,type="b",lwd=2,col="red")
This plot represents each of the sites in 2-D ordination space (x-axis = principal component 1, y-axis = principal component 2). (Should try to use something relevant for text)
Calculate and depict species loadings (i.e., principal weights in the eigenvectors) on each principal coordinate. Use the function envfit() along with the PC scores from our PCoA object. The function envfit() performs a linear correlation analysis based on standardized data (in other words, a simple linear regression) between each of the original descriptors (i.e., species) and the scores from each principal component. A permutation test is used to assess statistical significance, rather than using the F distribution.
print(vec.sp<-envfit(spe.pcoa$points, k=45, fert.data.4, perm=1000, na.rm=T))
##
## ***VECTORS
##
## Dim1 Dim2 r2 Pr(>r)
## Phylum -0.88963 -0.45669 0.8771 0.000999 ***
## Taxa 0.90418 0.42715 0.8888 0.000999 ***
## pH.experim 0.03946 -0.99922 0.2462 0.000999 ***
## Perc.Fertilization 0.22207 -0.97503 0.7852 0.000999 ***
## Insemination.mins -0.21279 0.97710 0.0960 0.002997 **
## Sperm.pre.exp.time -0.83941 0.54349 0.0583 0.041958 *
## egg.pre.exp.time -0.78719 0.61671 0.5652 0.000999 ***
## Sperm.per.mL -0.79476 0.60692 0.5732 0.000999 ***
## sperm.egg -0.99932 0.03686 0.6199 0.000999 ***
## n.females -0.66979 0.74255 0.5188 0.000999 ***
## n.males -0.68888 0.72487 0.2138 0.000999 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 1000
##
## 340 observations deleted due to missingness
p.max is the significance level that the species occurrence data must have with either PC in order to be depicted (these p-values were presented in vec.sp).
fert.data.4$pH.group <- cut(fert.data.4$pH.experim, c(6,7.3,7.5,7.8, 8.2))
pdf(file = "fert.PCoA.pdf", width = 9, height = 8)
pl <- ordiplot(spe.pcoa, type = "none", xlim = c(-1,1.5))
## species scores not available
points(pl, "sites", cex=0.8, pch=c(21,22,23,24)[fert.data.4$pH.group],
bg=c("red","blue","green","purple")[fert.data.4$Phylum])
plot(vec.sp, p.max=.01, col="black") # note, p.max = set p-value threshold for plotting vectors.
legend(x="topright", legend = levels(fert.data.3$Phylum), col=c("red","blue","green", "purple"), pch=c(16,16,16,16))
legend(x="right", legend = levels(fert.data.3$pH.group),pch=c(21,22,23,24))
dev.off()
## quartz_off_screen
## 2
dist.gower2 <- vegdist(fert.data.4[c(2,4:11)], "gower", na.rm = F)
spe.pcoa2 <- cmdscale(dist.gower2, k=2, eig=T, add=T)
print(vec.sp2<-envfit(spe.pcoa2$points, k=45, fert.data.4[c(2,4:11)], perm=1000, na.rm=T))
##
## ***VECTORS
##
## Dim1 Dim2 r2 Pr(>r)
## Taxa 0.63598 -0.77170 0.8539 0.000999 ***
## Perc.Fertilization 0.61708 0.78690 0.9546 0.000999 ***
## Insemination.mins -0.85462 -0.51925 0.0472 0.054945 .
## Sperm.pre.exp.time -0.98964 0.14356 0.0488 0.065934 .
## egg.pre.exp.time -0.99754 0.07008 0.5951 0.000999 ***
## Sperm.per.mL -0.99696 0.07787 0.6043 0.000999 ***
## sperm.egg -0.92696 0.37515 0.5829 0.000999 ***
## n.females -0.99770 -0.06773 0.5179 0.000999 ***
## n.males -0.99765 -0.06853 0.1923 0.000999 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 1000
##
## 340 observations deleted due to missingness
pdf(file = "fert.PCoA-nopH.pdf", width = 9, height = 8)
pl2 <- ordiplot(spe.pcoa2, type = "none", xlim = c(-1,1.5))
## species scores not available
points(pl2, "sites", cex=0.8, pch=c(21,22,23,24)[fert.data.4$pH.group],
bg=c("red","blue","green","purple")[fert.data.4$Phylum])
plot(vec.sp2, p.max=.01, col="black") # note, p.max = set p-value threshold for plotting vectors.
legend(x="topright", legend = levels(fert.data.3$Phylum), col=c("red","blue","green", "purple"), pch=c(16,16,16,16))
legend(x="right", legend = levels(fert.data.3$pH.group),pch=c(21,22,23,24))
dev.off()
## quartz_off_screen
## 2
All of the below except Sperm.pre.exp.time, sperm.egg, and n.males.
VECTORS Dim1 Dim2 r2 Pr(>r)
pH.experim -0.78537 0.61903 0.9157 0.000999 Perc.Fertilization -0.76812 -0.64030 0.8868 0.000999 Insemination.mins 0.42319 -0.90604 0.3439 0.000999 Fert.success.mins 0.38683 -0.92215 0.3583 0.000999 Sperm.pre.exp.time 0.95118 0.30863 0.0034 0.880120
egg.pre.exp.time 0.82583 0.56392 0.2503 0.001998 pH.delta -0.73572 0.67729 0.7674 0.000999 Sperm.per.mL 0.81881 0.57406 0.2492 0.001998 sperm.egg 0.20294 0.97919 0.0953 0.038961 *
n.females 0.90918 0.41640 0.2463 0.001998 ** n.males 0.86065 0.50920 0.0570 0.127872
pdf(file = "fert-correlation-panel.pdf", width = 12, height = 8.5)
pairs(na.omit(fert.data.4), lower.panel=panel.smooth, upper.panel=panel.cor)
dev.off()
## quartz_off_screen
## 2
Plot factors that were deemed sign. in PCoA
pH.experim -0.78537 0.61903 0.9157 0.000999 ***
# plot % fert ~ pH.experim by Phylum
ggplotly(fert.data.3 %>%
ggplot(mapping=aes(x=pH.experim, y=Perc.Fertilization, group=Phylum, col=Phylum, text=`Common name`)) +
geom_point(size=1.5, width=0.02) +
#facet_wrap(~Taxa) +
geom_smooth(method="lm", se=TRUE, aes(fill=Phylum)) +
ggtitle("Fertilization Rate ~ pH exposure by phylum"))
## Warning: Ignoring unknown parameters: width
## Warning: Removed 6 rows containing non-finite values (stat_smooth).
summary(aov(Perc.Fertilization ~ factor(Phylum)*pH.experim, fert.data.3)) #significant alone and interaction
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(Phylum) 3 15799 5266 7.590 5.82e-05 ***
## pH.experim 1 76312 76312 109.982 < 2e-16 ***
## factor(Phylum):pH.experim 3 10300 3433 4.948 0.00217 **
## Residuals 451 312932 694
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 6 observations deleted due to missingness
Insemination.mins 0.42319 -0.90604 0.3439 0.000999 ***
# plot % fert ~ pH.experim by Phylum
ggplotly(fert.data.3 %>%
ggplot(mapping=aes(x=Insemination.mins, y=Perc.Fertilization, group=Phylum, col=Phylum, text=`Common name`)) +
geom_point(size=1.5, width=0.02) +
#facet_wrap(~Taxa) +
geom_smooth(method="lm", se=TRUE, aes(fill=Phylum)) +
ggtitle("Fertilization Rate ~ Insemination minutes by phylum"))
## Warning: Ignoring unknown parameters: width
## Warning: Removed 65 rows containing non-finite values (stat_smooth).
summary(aov(Perc.Fertilization ~ factor(Phylum)*Insemination.mins, fert.data.3)) # not sign. alone, but sign for some phyla
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(Phylum) 2 22015 11007 14.131 1.18e-06 ***
## Insemination.mins 1 6 6 0.008 0.9288
## factor(Phylum):Insemination.mins 2 6310 3155 4.050 0.0182 *
## Residuals 394 306910 779
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 65 observations deleted due to missingness
egg.pre.exp.time 0.82583 0.56392 0.2503 0.001998 **
# plot % fert ~ pH.experim by Phylum
ggplotly(fert.data.3 %>%
ggplot(mapping=aes(x=log(egg.pre.exp.time), y=Perc.Fertilization, group=Phylum, col=Phylum, text=`Common name`)) +
geom_point(size=1.5, width=0.02) +
#facet_wrap(~Taxa) +
geom_smooth(method="lm", se=TRUE, aes(fill=Phylum)) +
ggtitle("Fertilization Rate ~ Egg pre-exposure minutes by phylum"))
## Warning: Ignoring unknown parameters: width
## Warning: Removed 266 rows containing non-finite values (stat_smooth).
fert.data.3$egg.pre.exp.time.log <- log(fert.data.3$egg.pre.exp.time+1)
summary(fert.data.3$egg.pre.exp.time.log)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 0.000 2.773 1.991 2.773 7.273 156
summary(aov(Perc.Fertilization ~ factor(Phylum)*egg.pre.exp.time.log, fert.data.3)) # not sign. alone, but sign for some phyla
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(Phylum) 2 2712 1356 1.911 0.149717
## egg.pre.exp.time.log 1 528 528 0.744 0.389184
## factor(Phylum):egg.pre.exp.time.log 1 8413 8413 11.855 0.000656
## Residuals 304 215740 710
##
## factor(Phylum)
## egg.pre.exp.time.log
## factor(Phylum):egg.pre.exp.time.log ***
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 156 observations deleted due to missingness
Sperm.per.mL 0.81881 0.57406 0.2492 0.001998 **
ggplotly(fert.data.3 %>%
ggplot(mapping=aes(x=log(Sperm.per.mL+1), y=Perc.Fertilization, group=Phylum, col=Phylum, text=`Common name`)) +
geom_point(size=1.5, width=0.02) +
#facet_wrap(~Taxa) +
geom_smooth(method="lm", se=TRUE, aes(fill=Phylum)) +
ggtitle("Fertilization Rate ~ sperm concentration (log-trans) by phylum"))
## Warning: Ignoring unknown parameters: width
## Warning: Removed 204 rows containing non-finite values (stat_smooth).
fert.data.3$Sperm.per.mL.log <- log(fert.data.3$Sperm.per.mL+1)
summary(aov(Perc.Fertilization ~ factor(Phylum)*Sperm.per.mL.log, fert.data.3)) # not sign. alone, but sign for some phyla
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(Phylum) 2 29032 14516 27.612 1.40e-11 ***
## Sperm.per.mL.log 1 537 537 1.022 0.313
## factor(Phylum):Sperm.per.mL.log 2 23492 11746 22.343 1.15e-09 ***
## Residuals 255 134059 526
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 204 observations deleted due to missingness
pH groups: – 6-7.3 – 7.3-7.5 – 7.5-7.8 – 7.8-8.2
ggplotly(fert.data.3 %>%
ggplot(mapping=aes(x=log(Sperm.per.mL+1), y=Perc.Fertilization, group=Phylum:pH.group, col=Phylum:pH.group, text=`Common name`, shape=pH.group)) +
geom_point(size=1.5, width=0.02) +
#facet_wrap(~Taxa) +
geom_smooth(method="lm", se=TRUE, aes(fill=Phylum:pH.group)) +
ggtitle("Fertilization Rate ~ sperm concentration (log-trans) by phylum"))
## Warning: Ignoring unknown parameters: width
## Warning: Removed 204 rows containing non-finite values (stat_smooth).
## Warning: Factor `Phylum:pH.group` contains implicit NA, consider using
## `forcats::fct_explicit_na`
## Warning: Factor `Phylum:pH.group` contains implicit NA, consider using
## `forcats::fct_explicit_na`
sperm.egg 0.20294 0.97919 0.0953 0.038961 *
# plot % fert ~ pH.experim by Phylum
ggplotly(fert.data.3 %>%
ggplot(mapping=aes(x=sperm.egg, y=Perc.Fertilization, group=Phylum, col=Phylum, text=`Common name`)) +
geom_point(size=1.5, width=0.02) +
#facet_wrap(~Taxa) +
geom_smooth(method="lm", se=TRUE, aes(fill=Phylum)) +
ggtitle("Fertilization Rate ~ sperm:egg ratio by phylum"))
## Warning: Ignoring unknown parameters: width
## Warning: Removed 2 rows containing non-finite values (stat_smooth).
summary(aov(Perc.Fertilization ~ factor(Phylum)*sperm.egg, fert.data.3)) # sign. main and interaction effects
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(Phylum) 2 16503 8252 29.06 1.16e-11 ***
## sperm.egg 35 135929 3884 13.68 < 2e-16 ***
## factor(Phylum):sperm.egg 1 9519 9519 33.52 3.09e-08 ***
## Residuals 180 51119 284
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 246 observations deleted due to missingness
n.females 0.90918 0.41640 0.2463 0.001998 **
# plot % fert ~ pH.experim by Phylum
ggplotly(fert.data.3 %>%
ggplot(mapping=aes(x=n.females, y=Perc.Fertilization, group=Phylum, col=Phylum, text=`Common name`)) +
geom_point(size=1.5, width=0.02) +
#facet_wrap(~Taxa) +
geom_smooth(method="lm", se=TRUE, aes(fill=Phylum)) +
ggtitle("Fertilization Rate ~ No. females used for assay, by phylum"))
## Warning: Ignoring unknown parameters: width
## Warning: Removed 97 rows containing non-finite values (stat_smooth).
summary(aov(Perc.Fertilization ~ factor(Phylum)*n.females, fert.data.3)) # sign. main effect, not interaction
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(Phylum) 2 14574 7287 9.605 8.62e-05 ***
## n.females 1 12866 12866 16.957 4.74e-05 ***
## factor(Phylum):n.females 2 2094 1047 1.380 0.253
## Residuals 362 274654 759
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 97 observations deleted due to missingness
summary(aov(Perc.Fertilization ~ factor(Phylum)*(pH.experim + Insemination.mins + egg.pre.exp.time.log + Sperm.per.mL.log + n.females), fert.data.3)) #
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(Phylum) 2 3959 1979 7.094 0.001069
## pH.experim 1 3416 3416 12.242 0.000582
## Insemination.mins 1 153 153 0.548 0.459852
## egg.pre.exp.time.log 1 1474 1474 5.284 0.022615
## Sperm.per.mL.log 1 4207 4207 15.079 0.000142
## n.females 1 7482 7482 26.815 5.68e-07
## factor(Phylum):pH.experim 2 995 497 1.782 0.171047
## factor(Phylum):Insemination.mins 1 10689 10689 38.308 3.64e-09
## factor(Phylum):egg.pre.exp.time.log 1 2883 2883 10.333 0.001536
## factor(Phylum):Sperm.per.mL.log 1 3397 3397 12.173 0.000603
## factor(Phylum):n.females 1 794 794 2.847 0.093164
## Residuals 190 53015 279
##
## factor(Phylum) **
## pH.experim ***
## Insemination.mins
## egg.pre.exp.time.log *
## Sperm.per.mL.log ***
## n.females ***
## factor(Phylum):pH.experim
## factor(Phylum):Insemination.mins ***
## factor(Phylum):egg.pre.exp.time.log **
## factor(Phylum):Sperm.per.mL.log ***
## factor(Phylum):n.females .
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 261 observations deleted due to missingness